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Abstract. Construction of splitting-step methods and properties of related non-negativity and 
boundary preserving numerical algorithms for solving stochastic differential equations (SDEs) of Ito- 
type are discussed. We present convergence proofs for a newly designed splitting-step algorithm and 
simulation studies for numerous numerical examples ranging from stochastic dynamics occurring in 
asset pricing theory in mathematical finance (SDEs of CIR and CEV models) to measure-valued 
diffusion and superBrownian motion (SPDEs) as met in biology and physics. 
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1. Introduction and examples. Stochastic differential equations (SDEs) are 
fundamental to describe and understand random phenomena in different areas of 
physics, engineering, economics, etc. Particularly, they serve as model for price fluc- 
tuations as in the famous Black-Scholes option pricing model, description of erratic 
movements of particles as in the Langevin equation or spatial processes like the super- 
Brownian motion. In most cases, explicit solutions of SDEs are very difficult to obtain 
and numerical approximations need to be exploited. In fact a wealth of methods for 
integrating numerically SDEs are known and tested ^| EH CHI 1271 QE| • 

In some practical applications autonomous Ito-type SDEs of the form 

(1.1) dX(t) = f(X(t))dt + a(X(t))dW(t), X{0)=X Q eV 

are not well defined unless a boundary condition is additionally given at the boundaries 
of the domain T> in which X{t) lives for all t (almost surely). For example, if X(t) 
is the price of a stock and l|l.l|) gives its time evolution, then V = [0, oo), where 
X(t) = implies an absorbing or reflecting state. Or, if X(t) models the number 
of certain species in a noisy environment, then T> — [0, if], where K is an attracting 
carrying capacity of the environment. 

As in most situations, boundary conditions are not needed to state the related 
problems (|1.1|) when the boundaries are unattainable in finite time. This is the case 
of natural boundaries, according to Feller's classification of one-dimensional diffusions 
[HI E3- The standard (unrestricted) Brownian motion on Hi is the most obvious 
example of diffusion with natural boundaries at infinity. The situation is different 
when the solution of i|l.l[) attains the boundaries in finite time. For example, the 
Brownian motion on IR + in which a boundary condition at x = needs to be specified 
to completely define the solution of it. Typical boundary conditions in this case are 
absorbing or reflecting ones and the solution of 1)1.1(1 depends on its specific choice, 
which is usually taken according to the nature of the problem under consideration. 

This problem of how to handle the boundary conditions also appears in the nu- 
merical approximations of (|l.l(l . where the naturally inherited boundary conditions 
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do need to be incorporated in the construction of numerical approximations. In fact, 
standard numerical methods such as Euler methods may fail to meet the boundary 
conditions, see |25j . This is also true for higher order converging Taylor-type methods 

(1-2) Y n+1 = ]T f a (Y n )I a , 

where a is a multiple index, TL a hierachical set of multiple indices, f a coefficient func- 
tions and I a iterated multiple stochastic integrals as derived from stochastic Taylor- 
type expansions (for its origin, see along time-discretizations 

= t < ti < t 2 < ... < t nT = T 

with step sizes A„t = t n +i — t n . For more details, see ^B]> 

Numerical time-discretizations of the type ((1.2(1 face two kinds of different (though 
related) problems when a boundary condition is specified. The first one is to restrict 
the values of the numerical approximations to live within the domain T>\ secondly, 
they also have to preserve the character of the boundary (natural, reflecting, ab- 
sorbing, etc.) in the numerical approximations of X(t). To exemplify the numerical 
approximation problems, let us consider the well-known square-root diffusion model 
of Cox-Ingersoll-Ross [3] 

(1.3) dX(t) = [a + bX(t)]dt + ay/X(t)dW(t), X(0) = x > 

with real parameters a > 0, b G IR and a > 0, which is widely used for interest rate 
modelling or as an alternative to geometric Brownian motion occurring in the Black- 
Scholes model of dynamic asset pricing in mathematical finance. It is well known 
that the (strong) solution of ((1.3|) is unique and preserves the non-negativity of the 
initial data. This property can be easily implemented in numerical approximations 
of ((1.3(1 by means of balanced implicit methods (see 12011221123); balanced Milstein 
methods ^B], composition methods based on Lie algebra techniques [21], but fail to 
incorporate the right boundary properties in the numerical paths. Other methods, 
while converging in the limit At — > 0, do not even conserve the non-negativity of the 
solution like the straightforward fixing by taking ^/|X(t)| instead of the last term in 
(Q|) (as done in [T2|1. 

The aim of this paper is to present an alternative strategy based on splitting- 
step methods to integrate numerically SDEs of the type Hl.lfl subject to boundary 
conditions, that both guarantee that the numerical solutions live in the domain T> and 
that the character of the boundary is preserved (a kind of numerically compact support 
property, i.e. the numerical approximation has to incorporate the fact that, for any 
positive X(t), there is a non-zero probability such that the stochastic process becomes 
zero at the next time-step). Moreover, the consistency of this new method shall be 
mathematically justified by two convergence theorems, namely one for sufficiently 
smooth Lipschitzian coefficients of involved SDEs fTheorem 13.1(1 and another, more 
complicated one for non-Lipschitzian SDEs (Theorem 13.2(1 . By those theorems we 
are going to establish the same convergence rates as standard, so far known methods 
possess. The second Theorem l3 . 2l covers the square-root case as exhibited by 1(1.3(1 too, 
and we present an alternative boundary- and positivity-preserving numerical splitting 
algorithm and its proof to that in ^21 without using substitutions such as y|X(t)| in 
diffusion terms of model 1(1.3(1 . 

The paper is structured as follows. After this introduction we introduce and 
discuss the fairly general splitting-step algorithm in Section 2. Section 3 reports on 
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general convergence theorems and its mathematical proof. Several numerical experi- 
ments in Section 4-6 strongly support the suggested splitting algorithm and results 
from previous sections by ordinary SDEs. In Section 4 we present some simple mod- 
els and related experiments. Section 5 is devoted to the simulation of processes as 
often met in dynamic asset pricing in mathematical finance. In Section 6 we men- 
tion some applications to measure-valued diffusions and reaction-diffusion equations 
demonstrating the potential range of the splitting-step algorithm to numerical treat- 
ment of random PDEs as well. Section 7 concludes with some supplemental remarks. 
Eventually, there is a small appendix on aspects of random number generation. 

2. Splitting-step algorithm. The general structure of a splitting-step algo- 
rithm, which is based on the idea in is as follows. Suppose that the more general 
equation which is to be integrated is of the form 

(2.1) dX(t) = [a(X(t), t) + f3(X(t), t)]dt + a(X(t), t)dW(t). 
We then decompose the above equation into the two equations 

(2.2) dXi{t) = P{Xx{t),t)dt + a{Xi{t),t)dW{t), 

(2.3) dX 2 (t) = a(X 2 (t),t)dt, 

where the splitting is done assuming that we know the exact strong solution for X\ (t) 
or the conditional probability V[Xi(t)\Xi(0)]. Thus, we can approximate the solution 
of i|2.1|l by a stochastic process Y t along time intervals [t, t + At] using the following 
two-step algorithm for each At, which we call splitting-step algorithm: 
Step 1. Knowing the value of Y t we obtain an intermediate value Y t which is obtained 

through the exact integration of (|2.2|l and Y t = X± (t + At) and with initial 

condition X\(t) — Yt. 
Step 2. Then Yt is used as the initial condition for (|2.3|l which is now integrated 

using any converging deterministic numerical algorithm to get X 2 . Then 

Y t+At =X 2 (t + At). 
The advantage of this splitting-step technique for SDEs subject to boundary 
conditions is is that if equation (|2.2|l is simple enough and we know the solution X\ 
of equation (|2.1J) . then the stochastic part of the problem can be handled correctly. 
For example, for the case (3(Xi,t) — and a(Xi,t) = \fX\, it is known that the 
conditional probability distribution is given by 

(2.4) V[ Xl (t)\xm] - \ (ffy) 1/2/ i (JV*i(*)*i(0)) e-*l*<*>+*<°» + 

+e-i x ^S(X 1 (t)), 

where I\ is the modified Bessel function and S(x) is the Dirac delta function. 

This cdf can be sampled using the rejection or inverse methods but this is compu- 
tationally expensive. However, we introduce here a very simple method |22j (see also 
U) for obtaining X\(f) by noting that the variable Z(t) = jX\{t) has a probability 
distribution given by the non-central ^-distribution, that is 

(2.5) V[Z\Z ] = ]T [ -^- r V xl . (Z) + e- x 'H{Z) 

i=i 3 ' 

where A = jZ and V x 2 (X) is the x 2 -pdf with 2j degrees of freedom. Equation l|2.5|) 
reveals that the probability distribution for Z is a linear combination of x 2 -pdfs with 
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Poisson weights. This fact can be exploited to generate X\ (t) efficiently. If we choose 
K from a Poisson distribution with mean A/2, then 



where Zi are independent Gaussian random numbers with zero mean and unit variance. 
Computationally, it is faster to sample the random number X)«=i z f using standard 
algorithms for random number generation of the x 2 distribution. Other examples of 
this sampling can be found in Section 0] and the Appendix. 

The obvious advantage of the splitting algorithm is that we exploit the structure 
of the original equation more efficiently. For example, in the above example, we get 
numbers from l|2.6Jl which are nonnegative and thus, if a(X, t) has nice properties, the 
approximated values for X(t) are nonnegative too. Such a splitting algorithm is not 
known from the literature to the best of our current knowledge, although the idea of 
splitting is not new for dynamical systems and their numerical integration. A different 
kind of splitting technique has been suggested in |13| . However, their algorithm called 
split-step Euler method is related to another subclass of splitting of SDEs and their 
resulting split-step algorithm is of lower order 0.5 of mean square convergence, which 
is restricted by their use of (drift-implicit) backward Euler method. Their method 
indirectly refers to the splitting 



of the original system in our set-up, where both equations for Xi and X 2 are 

numerically integrated in a separated fashion. In contrast to them, by a more effi- 
cient choice of splitting, we suggest even to remove stochastic integrals from numerical 
integration by appropriate splitting of f| 1 . l|l in order to achieve higher order of con- 
vergence and our technique is not only restricted to techniques of direct pathwise 
simulation. It may be noted that the explicit removal of stochastic integrals from nu- 
merical integration by splitting techniques is always possible and leads to converging 
numerical approximations with higher order under some appropriate assumptions on 
the diffusion coefficient a (such as uoi has sufficiently many bounded derivatives) in 
H 1 . Another form of splitting has been suggested by 0] . They present an algorithm 
which also takes advantage of techniques of numerical integration of ODEs. However, 
both 0] and |KSj do not discuss the issue of pathwise preservation of nonnegativity, 
monotonicity and boundedness by their numerical approximation techniques. 

3. General Convergence Theorems for Nonautonomous Equations. Let 

C* ,J (1R x [0, T]) denote the vector space of continuous functions / = f(x,t) which 
are i times continuously differentiable with respect to the space-coordinate Xk € M 
(k = 1,2,..., d) and j times continuously differentiable with respect to time-coordinate 



3.1. A first general convergence theorem. Recall that the original equation 

is 



(2.6) 




dXi(t) = [a(Xi(t),t) + 0(Xx(t),t)]dt, 
dX 2 (t) =a(X 2 {t),t)dW{t) 



te[o,T]. 
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Theorem 3.1. Assume that the coefficient functions a, € C 2,1 (JR d x [0, T}) and 
a G C 3 ' 2 (]R d x [0, T]) with exclusively uniformly bounded derivatives are such that 

sup IE [\X(t)\ 2 + \a(X(t),t)\ 2 + \(3(X{t),t)\ 2 + \a(X(t), t)\ 2 ] < +oo 

0<t<T 

for a fixed finite, nonrandom terminal time T > 0. Then the splitting-step algorithm 
with steps 1 and 2 has (global) strong and weak order 1.0 of convergence on the interval 
[0, T] (in the worst case). 

Proof. For simplicity, suppose that d = 1. Let 

= t < h < ... < t n < t n+l < ... < t N = T 

be any nonrandom partition of the given time-interval [0, T] with sufficiently small 
maximum step size 

A = max Iff — ti-\\ C 1. 

i=l,2,. .. ,N 

Define the local pathwise error by £^+1 = -^(^n + A„t) — X(t n + A n t) assuming that 
both exact solution X and its approximation X have started at the same value X(t n ) 
at time t n . To investigate this error, apply stochastic Taylor approximations to the 
processes X, Xi and X 2 - an idea which originates from the Wagner-Platen expansion 
|3"t] and was popularized by ^7] (more precisely speaking, this exploits the idea of 
an iterative application of Ito formula). For the sake of abbreviation, we shall write 
z n or z u for all occurring coefficients or processes (not referring to partial derivatives 
here), hence X n = X(t n ), a n = a(X(t n ),t n ), (3 n = f3(X(t n ),t n ), a n = a(X(t n ),t n ), 
similarly a u = a(X(u),u), f3 u = /3(X(u),u), o u = cr(X(u),u) and so on, unless it 
is stated differently wherever convenient. Furthermore, define the partial differential 
operators 

lVM = ^ + Mx,t) + nx,t)]V&) + ^(M)] 2 %^, 

Llf(x,t)=a(x,t)^p^, L\f(x,t) = Llf(x,t), 

where L° is mapping from C 2 » 1 (lR d X [0, T]) to C°<°(IR d x [0, T]) and L\ from C lfi (M d X 
[0,T]) to C°'°(]R d x [0,T]) for i = 1,2, and i!] from C 1 ' 1 ^ x [0,T]) to C°'°(lR d x 
[0, T]). First, apply stochastic Taylor approximations to the solutions of (|3.1|l to 
obtain 

X n+1 =X n + J " +1 [a(X(s),s) + p{X{s),s)]ds + J " + " a(X(s), s)dW(s) 
= X n + [a„ + [3 n ]A n t + a n A n W + 

rtn+i rs 

L° {a u + f3 u ) duds + I I L l {a u + (3 U ) dW '(u) ds + 



tn J t„ 
t n+ l 



L%a u dudW(s)+ / h\a u dW(u) dW(s) 



(3.4) = X n + [a n + (3 n ]A n t + a n A n W + ^Llcr n [{A n W) 2 - A n t] + i? ,„ 
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with remainder term 

rt n+ i fs 



°o,n — 



Lq{o: u + (3 U ) du ds + 

tn+l fS 

J Ll<j u dudW{s) - 



Ll(a u +0 u )dW(u)ds- 



tn + l fS PU 



LqLq<j v dvdW{u)dW{s) 



L L L L cr v dW(v)dW(u)dW(s) 

lt„ Jt n Jt n 

Now, we have to compare <|3.4[1 with what we get from the expansion of the splitting 
method. In the latter case, by application of Wagner-Platen expansion f3I] again, we 
arrive at 

X hn+1 =X n + f n+1 (3(X 1 (s), s )ds+ f n+1 a{X 1 (s),s)dW(s) 



X n + (3 n A n t + a n A n W ■ 

ftn+l fS 

L\j3 u du ds - 



tn 



tn + l rs 



t n +i rs 



L\f3 u dW{u) ds + 
r L\a u dW{u)dW{s) 



(3.5) 



L\a u dudW(s) 



X n + f3 n A n t + a n A n W + -L\a n [(A n W) 2 - A n t] + R u 



with the remainder term 

ftn+l r-s 

Ri,n = + / / L>i@ u du ds - 

Jt„ Jt n Jt n 

ftn+l r s 

L\a u dudW(s) 



tn+l 



L{/3 U dW{u) ds + 



tn+l fS pU 



L\L\a v dvdW{u)dW{s) 



t n + i PS pu 



L\L\o v dvdW(u)dW(s). 



Here, the coefficients (3 and a involved in the above integrals are evaluated at the 
arguments {X\(u\u) and {X\{v),v), respectively. Similarly, by deterministic Tay- 
lor expansion for the local approximation of (|3.3|) in the framework of the splitting 
method, one gets to 

/•tn+l 

(3.6) X 2 , n +i = Xi,n+i + J a(X 2 (s),s)ds = X 1<n+ i + a(Xi^ +1 ,t n )A n t + R 2 ^ n 



with remainder term 

ft n +l fS 



L^a u du ds 



tn+l fS - 



da{X 2 {u),u) da{X 2 {u),u) 
■ + a{X 2 (u),u)- 



du 



dx 



du ds. 



Here, the coefficients a involved in the above integrals are evaluated at the arguments 
(X 2 (u),u). An expansion of ot(Xi t n+i,tn) with respect to space- variable x gives 

(3.7) a(X hn+1 ,t n ) = a n + f + L^a(Xi(a), s) ds + f + L\a{X x {s),s) dW{s). 
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Now, plug expansions 1|3.7[1 and i|3.5[l into the expansion l|3.6|) in order to encounter 
with 

(3.8) X 2 ,„+i = X n + [a n + /3 n ]A n t + a n A n W + -L\a n [{A n W) 2 - A n t] + R 3 . n 



with the remainder term 

Rz,n — Ri.n + R2, 



LiOtsds 



L\a s dW(s) 



A„t. 



Consequently, the local pathwise error = X n+ i — X 2 , n +i can be represented by 



-n+l 



— X n ^ 



X2,n+1 — Ro,n — Rz,n 



tn + l ps 



t„ 



L^(a u + I3 U ) du ds 

tn+l l-s 

J L%a u dudW{s)-\ 

tn+l ps pu 



tn + l pS 



t„ 



Ll(a u + (3 U ) dW(u) ds + 



tn 



tn+l ps pU 



LlL\a v dvdW(u) dW(s) 



LlL\a v dW{v) dW(u) dW{s) 



Llp u duds - I I L\f3 u dW(u) ds 

t„ Jt n Jt n 

tn + l p s ptn+1 P s P u 

LP x a u dudW{s)- 



L\L\o v dvdW(u) dW(s) + 



tn + l pS pU 



L\L\a v dvdW(u)dW(s) 



tn + 1 pS 



L%ol u du ds + 



tn 



(3.9) 



L\a s ds 



L\a s dW(s) 



A n t. 



Now recall that a, /?, a have exclusively uniformly bounded derivatives and their sec- 
ond moments along the solution of (|3.1(l are bounded on [0, T] . Therefore, all operators 
L\ applied to coefficients a, (3 and a have images with uniformly bounded second mo- 
ments. This implies that = 0{[A n t} 2 ) and E {e l ° c - TEe l ° c ) 2 = 0{[A n tf). For 
mean square convergence, it remains to apply fairly general convergence theorems as 
known from or with local rates r\ — 2 and r 2 = 1.5 in order to conclude 
the global strong (i.e. in the £ 2 -sense) convergence order r = 1.0. More precisely, we 
have 



max [JE\X(t n )-X(t n )\- 

=0,1,. ...N y 



1/2 



= 0(A) 



provided that the initial errors IE \X(0) - X(0)\ 2 = 0{A 2 ) started at L 2 -integrable 
initial values which are independent of the cr-algebra generated by the underlying 
Wiener process W . Eventually, the global order r w = 1.0 of weak convergence with 
respect to smooth, polynomially bounded test functions is rather obvious from the 
fact that the local weak rate r w — 2.0 (exploiting standard techniques known from jl()| 
and 53 on weak convergence analysis). Consequently, the splitting-step algorithm 
has the claimed convergence orders under the above stated assumptions, and the proof 
is complete. □ 
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3.2. Numerical generation of strong path solutions for some processes. 

The success of the splitting-step method requires the exact numerical integration of 
step 1 which can be problematic in some specific situations. That is why we have ap- 
plied in previous numerical examples such equations which allow us to represent the 
solution from step 1 in terms of pure functions of the Wiener process and initial values 
such as X(t) — F(X(Q),W(t)) or simple deterministic integrals as seen with the ex- 
ample of Bessel-type diffusions X(t) = (y/X(0) + W(t)) 2 or the geometric Brownian 
motion X{t) = A(0)exp((a - a 2 /2)t + aW(t)). Also, for the Ornstein-Uhlenbeck 
processes or logistic equations, one may use well-known integration-by-parts formula 
and / or the information on the exact distribution of involved stochastic integrals with 
deterministic differentials which can be pathwisely treated by deterministic quadra- 
ture methods. As another alternative, one could exploit the Doss representation (see 
[S]) to develop a semi-analytic approach or a ODE-PDE approach to decompose the 
original problem into a random ODE and a deterministic PDE in order to figure out 
which numerical implementation is more efficient in conjunction with our splitting 
technique. However, in general one produces additional discretization errors which 
might influence the accuracy of the computations using the splitting algorithm as 
well. Then it is a must to consider its stability properties too. Such a fairly complex 
and very problem-dependent work we leave to the future. 

3.3. A general theorem on L 2 -convergence based on VOP. To relax some 
of the technical assumptions, once may also exploit the variation-of-constants formula 
(VOP). Suppose that the original equation is 

(3.10) dX(t) = [a[X[t),t) + (3(X(t),t)]dt + a 1 (X(t),t)dW 1 (t) + <r 2 (X(t),t))dW 2 (t) 
where W\ and W2 are independent Wiener processes. Consider the splitting 

(3.11) dXi{t) = + ai(X 1 (t),t)dWi(t), 

(3.12) dX 2 {t) = a(X 2 (t),t)dt + a 2 (X 2 {t),t)dW 2 (t). 

Let Ci ocLip {S) denote the vector space of local Lipschitz continuous functions on the 
open set S. 

Theorem 3.2. Assume that the coefficient functions a,fie Ci ocLlp (ID x [0, T\) 
and o~i £ Ci ocLip (IDx [0,T]) are such that the continuous unique strong solution X of 
\3.1U\) exists on the closed set ID, 



sup IE 

0<t<T 



\X(t)\ 2 + sup sup IE |$(s,A(i))| 5 

- 1 0<t<Tt<s<T L 



< 



sup IE 

0<t<T 



\[^(s, X(t))]~ 1 a(X(s), s)\ 2 ds 



[$( S ,A(t))]-V 2 (A( s ), s )| 2 d s 



<+oo 



for a fixed finite, nonrandom terminal time T > and the stochastic flow $ generated 
by is mean square Holder-continuous with exponent rn > 0.5. Furthermore, 

suppose that step 1 of the splitting algorithm is exactly integrable and step 2 can be 
carried out with local mean accuracy with rate r% > 1.0, local mean square accuracy 
with rate r 2 > 0.5 and 



miniji, rn + 1.0) > min(r 2 , rn + 1-0) + 0.5. 
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Then, in the worst case, the error of splitting-step algorithm with steps 1 and 2 has 
(global) order 

r g > mm(1.0, min(r 2 , ru + 1-0) — 0.5) 

of L 2 -convergence on the interval [0, T]. 

Proof. Suppose that (|3 . 1 II) has known fundamental solution $ = $(t, Xq) with 
a.s. Holder exponent th- Then the exact solution of the original equation (|3.10(l 
possesses the pathwise representation 

/t+At 
[$(s,X(t))]" 1 a(X(s),s) ds 



-$(t + At,X(t)) 



t+At 



[$(s,X(t))]-V 2 (X(s),s) dW 2 (s) 



on each subintervals [t, t + At] C [0,T]. Note that $(t + At,X(t)) is independent 
of W 2 (s) — W 2 (t) for s > t. Let Y(t) denote the value of right-continuous numeri- 
cal approximation of splitting step algorithm (which we always can construct using 
step functions in a standard way). The main idea is to apply the fairly general L 2 - 
convergence theory known from ED and |2H]. For this purpose, we need to 

study the local conditional accuracy of our splitting algorithm. Locally, we may sup- 
pose that X(t) = Y{t) = x is deterministic (jFt-adapted). First, consider the local 
conditional mean accuracy of the splitting algorithm. We find that 



IE 



X(t + At) - Y(t + At) 

t+At 



Tt 



IE 



$(t + At, x) 



[<P{s,x)]- 1 a(X(s),s)ds+ 



< 



-<f>{t + At,x) 

rt+At 



l + AI 



IE 



IE 



IE 



t+At 



[$(s,x)}- 1 a 2 (X(s), s) dW 2 (s) - {Y{t + At) - x) 
$(t + At, x) - $(s, xfj [$(s, x)r 1 a(X(s), s) ds 
($(t + At, x) - x)\ x)]- 1 a 2 {X{s), s) dW 2 (s) 



T 



T 



T 



t+At 



t+At 



T 



a(X(s),s)ds + J a 2 (X(s),s)dW 2 (s)-(Y(t + At)-x) 

< ^ 1 [At] mm ( ri '( 2r "+ 1 )/ 2+0 ' 5 ), 

where K\ is a real constant, hence the local conditional mean accuracy rate 

min(r\, (2r# + l)/2 + 0.5) = min{r\, ru + 1-0) 

of the splitting algorithm can be verified. Second, consider the local conditional mean 
square accuracy of the splitting algorithm. 



IE 



|X(t + At) -F(t + At)p 



T 



1/2 
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= IE 



$(t + At,x) 



t+At 



[$(s,x)]~ 1 a(X(s),s) ds+ 



t+At 



1/2 



+$(< + At, x) / [*(s, x)]- 1 a 2 (X(s), s) dW 2 (s) - (y(i + At) - x) 



< IE 



IE 



t+At 



t+At 



($(t + At, i) - $(s, a;)) [$(s, x)]" 1 a(X(s), s) ds 



($(t + At, x) - $(s, au)) [$(«, x)]- 1 a 2 (X(s), s) dW 2 (s) 



1/2 



1/2 



IE 



t+At 



t+At 



1/2 



Tt 



a(X(s),s)ds + J a2(X(s),s)dW2(s)-(Y(t + At)-x) 

< K 2 [At] min( ' r2 - ( - 2rH+1 ' )/2+0 - 5] , 

where if 2 is a real constant, hence the local conditional mean square accuracy rate 

mm(r 2 , (2m + l)/2 + 0.5) = min{r-i,rn + 1.0) 

of the splitting algorithm can be derived. Now, apply the general convergence theory 
from JU], [27| and [2H| to conclude the worst case estimate of global mean square 
convergence rate r g as stated in Theorem 13.21 along nonrandom partitions of [0, T]. 
This completes the proof of Theorem 13.21 □ 



Remark. Consider the example 

(3.13) dX{t) = (1 + X{t))dt + 2y/X(t)dW(t) 

with a(X,t) = X, (3(X,t) = 1 and ai(X,t) = 2\/X and a 2 (X,t) = 0. This example 
satisfies the assumptions of Theorem 13.21 with rn = 0.5. To see this fact, one has to 
show that the related stochastic flow is given by the Bessel-type flow 

<f(t,X(0)) = (^X{0) + W(t)) 2 

while using Ito formula. Moreover, this flow has uniformly bounded moments of 
all orders for nonrandom initial data and is mean square Holder-continuous with 
exponent r# = 0.5 since 

$(t, X(0)) - $( s , X(0)) = 2y/X@j(W{t) - W(s)) 

and 

IE |$(t, X(0)) - $(s, X{0))\ 2 = 2IE [X(0)]IE \W(t) - W(s)\ 2 = HE [X(0)}\t - s\ 

provided that X(Q) > (a.s.) is independent of the process W. Therefore, the related 
splitting-step algorithm can achieve a global L 2 -convergence order r g = 1.0 since the 
numerical integration of step 2 can be implemented by any deterministic Runge-Kutta 
method with an interplay of local accuracy rates r\ — 2.0 and r 2 = 1.0. Similarly, 
we can proceed for other equations with \fi^)- or other Holder-continuous terms with 
Holder exponent > 0.5. 



NON-NEGATIVITY PRESERVING NUM. ALGORITHMS FOR SDEs 



11 



>- 

O 10 r 
C 



10 - 



p 



10" 



10 
At 



10" 



Fig. 4.1. Absolute value of the error \4-Hj as a function of At for t = 1.0 and X(0) 
obtained using the splitting- step algorithm. The dashed line is proportional to At. 



4. First illustrative numerical experiments. In this section we will give 
several illustrative and important examples of applying our new stochastic numerical 
schemes to SDEs and test the rate of convergence obtained in previous Section [3] 

4.1. Using the transition probability. We applied the splitting-step method 
to the following SDE 



(4.1) dX(t) = (1 + X(t))dt + 2 v / X(t)dW(t) 

with a(X,t) = 1 + X, P(X,t) = and cr(X, t) — 2\/X. The conditional mean value 
is given by 

(4.2) IE(X(t)|X(0) = a;o) = (2;o + l)e t -l 
for any non-random value xq > 0. In figure |4~T1 the error 

(4.3) si = |]EA > (i)-[(xo)e t -l]| 



versus decreasing uniform step size At is depicted, where X(t) is the solution obtained 
through the numerical approximation. Figure l4~T1 shows statistical-numerical evidence 
that the method has weak order 1.0 while using constant step sizes At. 

4.2. Using the exact solution. We now apply the splitting-step algorithm to 
the stochastic Ginzburg-Landau equation 

(4.4) dX(t) = (X(t) - [X{t)} 3 ) dt + X(t) dW(t). 

In this case we take advantage of the known exact solution for the linear part of this 
equation which is 

(4.5) dX x {t) = Xi{t)dt+X x {t)dW{t) => Xy{t) = Xi(s)exp ((t-s)/2+W(t)-W(s)). 

To integrate the remaining part of the equation we use an partial-implicit nonstandard 
technique which is a nonnegativity-preserving one 1 given by 

(4.6) X 2 (t + At) = X 2 (t) - ^Xi(t)[X 2 (t) +X 2 (t + At)] 

1 This is true only if [X2(0)] sup ngJV A„t < 2. 
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At 

Fig. 4.2. Value of the error \4- 7| ) as a function of At for t = 5.0 and X(0) = 1.0 using the 
splitting-step and Euler algorithms. Dashed lines are the power laws At and V At. 

Our results for the strong error 

(4.7) e k (t) = {M\X(t)-X(t)\ k ) 1 / k 

are shown in Figure l4~2l for k — 1,2 and compared with the Euler algorithm for the 
same equation. Obviously, our results indicate that our proposed splitting method 
indeed is of strong order 1.0 while using constant step sizes At. 

5. Stochastic models in finance. In this section we apply the splitting-step 
method to some fundamental models in mathematical finance. This application will 
also serve to introduce other algorithmically simple samplings of conditional proba- 
bility transitions as for the \J X\(t) case shown in the introduction. 

5.1. Interest rate model of Cox-Ingersoll-Ross. An interesting example in 
which the splitting-step scheme is particularly efficient is the Cox-Ingersoll-Ross (CIR) 
model [3] for stochastic interest rates 

(5.1) dX(t) = [a+ bX(t)]dt + ay/X(t)dW(t), X(0) =x >0 

with real parameters a > 0, b G IR and a > 0. As stated in the introduction, 
strong solutions of (|5.1|) are nonnegative. However, depending on the specific values 
of parameters a, b and er, distinct behavior at the boundary is possible. 

• If a > a 2 1 2 then the solution is always positive X(t) > if x > 0, because 
the boundary X(t) — becomes unattainable. 

• If a < a 2 / 2 there are infinite many values of t > for which X(t) = 0. The 
boundary becomes attainable, but it is (instantaneously) reflecting. That is, 
when a sample path reaches 0, then it returns immediately to the interior of 
the state space in a reflecting manner. 

The exact transition density for the CIR process is known, but its sampling can be 
difficult depending on the parameters a, b and a. Here we may exploit the simplest 
splitting which appropriately reflects the boundary behavior of the CIR process by 
SDEs 




(5.2) 
(5.3) 



dXi(t) =adt + ay/xJf)dW{t) 
dX 2 (t) = bX 2 (t)dt 
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Fig. 5.1. Numerical approximation of the CIR process 15. 1( 1 using the splitting- step method 
15. SK . (solid line) and the method proposed in (dashed line) with a = b = l and At = 10 -2 and 
X(0) = 1 and different values of a: Upper panel shows the case a > <r 2 /2 where the boundary at zero 
becomes unattainable, while lower panel is for a < cr 2 /2. The splitting-step numerical approximation 
reproduces the attainability of the boundary condition while preserving the non-negativity of the 
solution of 15. 1H . whereas the method in MSSI has trajectories which take negative values for X(t). 



which can be easily inferred by noting that the boundary behavior does not depend on 
the parameter b. To integrate the first step in the splitting-step system H5.2|) - (j5.3(l we 
note that the process defined by 15.2|l represents an a-dimensional Bessel process . 
Its transition density V \X\ (t + At) \Xi (t)] can be written in terms of a non-central \ 2 
distribution and in particular, we have that: 

(5.4) Xl {t + M) = ^^{\) 



where 

(5.5) A 



4Xi (t) 4a 
a 2 &t ' a 2 



and x'dW random numbers can be sampled using the algorithms in the appendix. 
The last part of the splitting-step scheme l|5.3|l can be integrated exactly or using 
numerical approximations. In our case we use deterministic Euler approximations 
where non-negativity is preserved if At is small enough. Our simulations for the CIR 
process are shown in Figure l5~Tl The boundary behavior for different values of a and 
a is reproduced and, at the same time, non-negativity is conserved. Note that other 
usual integration strategies as that of |12| eventually produce negative values for X(t) 
(for finite At) which lack any possible interpretation in the context of finance and 
could induce severe errors in option valuation. 

5.2. Constant Elasticity Volatility models. Another important stochastic 
process in finance is the constant elasticity of variance (CEV) diffusion to model asset 
prices. This process, first introduced to finance by Cox 0, is capable of reproducing 
the volatility smile observed in the empirical data unlike other standard price models 
like the Black-Scholes-Merton geometric Brownian motion. The process is defined as 



(5.6) dX(t) = nX(t)dt + aX(tydW(t), X(0) = x . 
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The CEV model includes the geometric Brownian model of Black, Scholes and Merton 
(7 = 1) and the square-root models of Cox and Ross (7 = 1/2). Contrary to the Black, 
Scholes and Merton model, the CEV model incorporates a variance adjustment that 
causes the absolute level of the variance to decline as the stock price rises and to 
rise as the stock price declines, as seen empirically in most equity and interest rate 
volatilities. 

For our purposes, the CEV model encompasses most of the boundary conditions 
we can implement at X(0) — 0. Depending on the value of 7 the boundary X(0) = 
is 

• Natural boundary for 7 > 1, which means that boundary is unattainable in 
finite time. 

• Exit or absorbing boundary for 1/2 < 7 < 1, i.e. X(t) reaches zero with finite 
probability in finite time and gets absorbed at it. 

• Regular boundary for 7 < 1/2. Now the boundary can be reached in finite 
time, and we need to specify a boundary condition at X(0) = 0. Typical 
choices are reflecting or absorbing ones. 

The transition density for the CEV process is known for general values of 7 and fj, (see 
0). However, due to the fact that the boundary behavior of X(t) does not depend 
on /1, we propose here the following splitting-step system 

(5.7) dX x = aXjdW{t), 

(5.8) dX 2 = fiX 2 dt. 

The classification of the boundary Xi(t) = for (|5.7|l is the same as for l|5.6|l . We 
discuss now the sampling of the transition density for each value of 7. 2 

5.2.1. 7 > 1, Natural Boundary. In this case X\(t) — is unattainable and 
its probability transition density is given by 



(5.9) V[x t \x ] 



1/2-27 1/2 
H x o 

^(7-1) 



■ exp 



2(1-7) 
^0 



+ X 



2(1-7) 



2a 2 (l- 7 ) 2 £ 



I 



1— 7 1 — 7 



a 2 (l- 7 ) 2 i 



To sample this distribution, we can consider the relationship between the CEV and 
the CIR processes. Since Xi(t) = is not accessible, we can make the following 

,2(1-7) 



change of variables Y x = X^ 1 ' /[4:(j - l) 2 cr 2 ] to find that dY x = Xdt + y/¥[dW(t) 
where A = (1 — 27)/[4(l — 7)]. Thus Y\ is a Bessel process which can be sampled 
using the non-central x 2 distribution. In the end we have that 



(5.10) 



Xx(t + At) = 



(7-l)VAt^ 2 



Xi(t) 



2(l- 7 ) 



o- 2 (7- l) 2 At 



where d = (1 — 27)/(l — 7) > 1/2. A typical path of the CEV process for 7 > 1 is 
shown in Figure l5~2*l 

5.2.2. 1/2 < 7 < 1, Absorbing boundary. In this case X\(t) = is an 
absorbing boundary. The transition probability density was found by Cox and is 
given by 



(5.ll)V[x t \x ] 



x 



1/2-27 1/2 



a 2 (l- 7 ) 



exp 



x, 



2(1-7) 




+ X 



2(1-7) 



2a 2 (l- 7 ) 2 t 



^1 — 7 1— 7 



CT 2 (l- 7 ) 2 t 



J We do not consider the case 7 = 1, since this is trivially integrated using the strong solution 
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Fig. 5.2. Numerical approximation of the CEV process 15. 61) using the splitting- step method 
j5.7f - f5.HII . for different values of 7. The case 7 = —1 is the solution of 15, 61) with reflecting 
boundary conditions, while the case 7 = 3/4 gets absorbed at zero at finite time (not shown in the 
logarithmic scale). Parameters are p. = 0.1, At = 10~ 2 . 



This transition probability does not integrate to one because there is a finite proba- 
bility that the trajectory gets absorbed at zero given by [2] 



(5.12) 



VMXtit)] = G 



Xl{t) 



-2(1-7) 



2(1-7)' 2cr 2 (l-j) 2 At 



where G(v, x) is the complementary Gamma function. For this equation, we know that 
trajectories of (|5.7|l get absorbed at zero with probability one while taking the limit 
At — > 00. Actually, this is also the case for the solutions of l|5.6|l for 1/2 < 7 < 1. This 
is counterintuitive since the mean IE [Xi(<)] is constant for 1)5. 7f) or grows exponentially 
like x^e^ 1 for 1)5. This intriguing feature of l|5.t)[) hampers the numerical simulations 
of this process and in fact only exact sampling of the probability transition density 
H5.11fl and (|5.12|) correctly accounts for it at finite At. 

While similar to l|5.9[) . the transition probability density (|5.11() cannot be sampled 
using the relationship to the CIR process, since the solution of dX\ = XjdW(t) has 
a finite probability to be absorbed at Xi(t) — 0. However, using the relationship 
I n (x) — I- n (x), n € N for the modified Bessel function we have that if 1/[2(1 — 7)] = n, 
X\ can be sampled from a non-central % 2 distribution with negative (integer) number 
of degrees of freedom (see appendix) : 



(5.13) 



Xt(t + At) = 



(7-l)VAi x « 



2(1-7) 



a 2 (7- I) 2 At 



where d = 2 - 2n = 0,-2,-4, and then 7=1- l/2n = 1/2,3/4,5/6, .... In 

figure IS~2l we show a typical path for the CEV process for 7 = 3/4 (d = —2) in which 
we see that it gets absorbed at zero for finite time. 

5.2.3. < 7 < 1/2, Regular Boundary. Now the boundary is attainable 
and, contrary to the previous case, a boundary condition should be provided. If we 
consider the simplest cases, namely absorbing or reflecting boundary condition then 
the transition probability density is given by 



{5.l4?±[x t \x ] = 



1/2-27 1/2 
0*11 - 7 | 



■ exp 



2(1-7) , 2(1-7) 
b _0 ~r x t 

2(7 2 (1 -j) 2 t 



1— 7 1— 7 
x x t 

o 2 (l - 7 ) 2 t 
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where V+ corresponds to the absorbing case and V- to the reflecting one. In the 
latter case, the relationship of the CEV and CIR processes can be exploited to sample 
the probability distribution Q5.14[) . with the same result as in (|5.1U|) . Simulation of 
the CEV process in this case is shown in figure (|5.2|l . However, in the absorbing 
case, V- cannot be sampled either by using the corresponding CIR process or by 
using negative dimensions as in the previous case, since for 7 < 1/2 we have that 
1 < l/[2(7 — 1)] < 0. Thus V- should be sampled using rejection or transformation 
methods. 

6. Measure valued diffusions / Reaction-diffusion problems. During the 
last decades much attention has been devoted to stochastic spatial models of interact- 
ing and branching particles like the contact model, the voter model, the normal and 
oriented percolation, etc. 1301 EJ 03 Despite its simplicity these models display inter- 
esting critical properties at high spatial dimensions and serve as universality classes 
for more complicated situations. While much analytical progress has been reached 
in the study of this models, some properties of them must be understood by using 
numerical methods. In this respect, several efficient particle stochastic simulations 
has been proposed and studied. An alternative approach is to study the convergence 
of the particle process to a continuum measure value diffusion in which the system 
is described by the (stochastic) concentration of particles p(x, t) at a given spatial 
location x. This accelerates the numerical simulations and also could help to iden- 
tify the relevant dynamics at the coarse-grained dynamics. Useful representations 
of these measure-valued diffusion are the ones interpreted as solutions of a martin- 
gale problem in terms of stochastic partial differential equations. Despite its clear 
interpretation, this representation is not usually considered in numerical simulations. 
The reason for that is the inability of usual numerical methods to handle correctly 
the non-negativity character and the Poissonian fluctuations of the concentration of 
particles close to p — 1221^3- ^ n this respect, we will see that the splitting scheme 
provides a very efficient and accurate method to study these models. 

6.1. Super-Brownian motion. The most simple and studied measured- value 
diffusion is the super-Brownian motion. The super-Brownian motion arises as the 
scaling limit in various critical branching systems when the interaction between them 
is weak, i.e. when the system is above some critical spatial dimension . Above this 
critical dimension we expect a Gaussian limit and indeed the super-Brownian motion 
is the Gaussian limit of a number of models: the voter model above 2 dimensions, 
the contact process above 4 dimensions, oriented percolation above 4 dimensions and 
percolation over 6 dimensions. 

Super-Brownian motion can be studied analytically by using the log-Laplace 
transform that maps its dynamics into a non-linear partial differential equation, a 
result due to Dynkin [Hj. More general situations or particular properties of the sBm 
can only be reached through numerical simulations. To our knowledge there is no 
numerical simulation of the martingale problem of the sBm. In one dimension, the 
martingale problem of the super-Brownian (sBm) motion is described by the stochas- 
tic partial differential equation 

(6.1) dp{x, t) = Ap(x, t)dt + y/ap(x, t)dW(x, t) 

where W(x, t) is a Wiener sheet. It is well known that the solutions of the sBm 
die in finite time almost surely. Another interesting property is that the support of 
the solution, i.e. the set for which u[x, t) > is compact, provided that the initial 
condition has a compact support. 
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Fig. 6.1. Realization of the sBm in one dimension. The figure shows density plots for u(x,t) 
as a function of time with initial condition u(x, 0) = 0.1. The solid line depicts the extremes of the 
support at each time. Parameters are Ax = 1, At = 0.1 and o = 1.0 in a L = 128 lattice. 



We use the splitting-step method for approximating the strong path solutions of 
1)6.1(1 . To this end we approximate the sBm by the super-random walk on Z d (see [H]) 



(6.2) dUi{t) = A iUi (t)dt+ ^^0±dWi(t) 

where A^ is the discrete Laplacian operator 

u i+1 - 2ui + 



Am = 



(Ax) 



in Z d , Ax is the lattice spacing and Wi(t) are independent Wiener processes in time t. 
Equation ((6.2(1 represents a model of interacting Feller diffusions. The splitting-step 
algorithm in this case is based on 

(6.3) d4\t) = ^0W i(t ) 

(6.4) duf\t) = A i uf ) 

where the first step is integrated using the transition probability 1(2.5(1 and its sampling 
1(2.6(1 and the equation below can be integrated using standard schemes. In Figure 
1(6.1(1 we observe a simulation of the sBm in one dimension. Our method not only 
provides an accurate description both in the strong and weak sense of the sBm, but it 
also incorporates one of the main properties of the sBm, namely the compact support 
property. 

In two dimensions equation ((6.1(1 is not well defined, but still we can study equa- 
tion 1(6.2(1 in the lattice 1? . As our simulations show, the compact support properties 
of sBm in two dimensions are preserved and we also see the cluster formation and 
their disappearance at large times. In contrast to that, the methods of Gaines do 
not provide such efficient numerical approximations. 
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Fig. 6.2. Strong path approximation of the sBm in two dimensions. The figures show density 
plots of u(x,t) surrounded by the border of the support (solid line) at different times. The initial 
condition is u(f,0) = 0.1. Parameters are Ax = 1, At = 0.1 in a 256 X 256 lattice. 



6.2. Contact process. The contact process is a model of spreading of infection 
in a lattice in which an site can be infected via contact with an infected site in its 
neighborhood 30 . Infection and recovery happens at different rates and there is a 
critical value of them for which an infection started from a single infected individual 
will die out in finite time. The contact process with finite range of infection has a 
critical dimension d c = 4. At higher dimensions the critical contact process converges 
to the sBm, but not below d c . At lower dimensions, it is not known what the scaling 
limits should be. 

However, long-range interaction with suitable scalings show that the contact pro- 
cess converges to the sBm for d > 2, see [7]. In one dimension, Mueller and Tribe 
showed that the rescaled density of particles for long-ranged contact process weakly 
converges to the solution of the SPDE 

(6.5) dp(x, t) = [Ap(x, t) + 0p(x, t) - p(x, tf]dt + y/pdW(x, t) 

where W(x, t) denotes a space-time Wiener process. Moreover, they showed that the 
above equation undergoes a phase transition at a critical value of 9 C for which 



(6.6) V(u(x,i) survives) 



= if 6 < 
> if 6 > 



The non-trivial behavior of the solution of l|6.5|l is believed to represent the well known 
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Fig. 6.3. (left) Strong path approximation of equation 1 6. 51) for Ax = 1 and L = 2 14 and 
different values of 9 below and above the phase transition. (Right) Critical value of 9 C as defined in 
16.61 ) as a function of At. The straight line is a linear fit to the data. Parameters are Ax = 1,L = 
2 14 . 



universality class for contact processes (also named directed percolation universality 
class 0E21). 

As before, numerical simulations of (|6.5|) can be now addressed using the splitting 
scheme and the efficient random number generators for the conditional probability 
|22|. In particular, we discretize the spatial operators in a lattice Z [like in JOJ] and 
split the dynamics as follows: 

(6.7) dpf } = sfpY>dW{t) 

(6-8) dp^ = [A iP f^ep^-{pf^]dt 

where the last equation is numerically integrated using Euler approximations with 
sufficiently small step sizes. Results for strong approximations of pi(t) are shown in 
Figure 16.31 where we can see a typical realization for the subcritical (infection dies 
out), critical and super-critical (infection spreads) for different values of 9. 

We can calculate the critical value of 9 C in one dimension using our algorithm. 
To this end, we identify the critical point 9 C using equation 1)6. 6f) and finite scaling 
techniques of statistical mechanics. In particular we found that 9 C = 0.777 ± 0.001 for 
Ax = 1, and the convergence to this value is of order one. 

7. Conclusion. The general idea of this paper is to propose a splitting of SDE 
Ijl.lfl into two new SDEs for which one can keep nonnegativity during integration 
of both subsystems. This is achieved by either solving one subsystem exactly in 
the pathwise sense or using its transition probabilities, and solving the other (here 
nonrandom) subsystem by nonnegativity-preserving numerical methods. In this way 
one is able to preserve nonnegativity and a maximum of convergence order 1.0 both 
in weak and strong sense. For the efficiency of our splitting algorithm, it is crucial 
to find a splitting into appropriate subsystems. For this purpose, one also tries to 
incorporate more complicated boundary conditions of the original SDE p. 1(1 into an 
explicitly solvable subsystem such that a fairly easier numerical integration of the 
remaining subsystem remains to be implemented. For example, for SDEs 



dX(t) = [a(X[t),t) + XX(t)]dt + aX(t)dW(t), 
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one makes use of the splitting into 

dXi(t) = \X x (t)dt + aX x (t)dW(t), 
dX 2 (t) = a(X 2 (t),t)dt, 

where the explicit solution of the first component X\ is given by 

X x {t) = Xi(0) exp((A - X -G 2 )t + aW(t)) 

which possesses the monotone property of leaving the positive axes [0, +00) invariant 
by this type of random mapping (almost surely). This idea can be easily extended to 
nonlinear systems of SDEs with its splitting into linear and nonlinear subsystems in 
several dimensions. Another type of splitting is found for nonlinear systems 

(7.1) dX{t) = f(X(t),t)dt + a(X(t), t)dW(t) 

with continuously differentiable coefficients a as follows. Rewrite this equation to as 



dX(t) 



dt 



2 v w ' ' dx 2 v w ' ' dx 

+ (T{X{t),t)dW{t) 

and set 

a{x, t) = f(x, t) - ^(r(x, t) ^ , 

nl , 1 , .da(x.t) 
P(x,t) - -a(x,t) ^ > 

Then the splitting-step algorithm is applied to system (JTi,^) satisfying equations 
(12. 2|) and (|2.3|) with coefficients a and (3 as defined above. This works at least ef- 
ficiently if cr{x, t) = a(x) does not depend on t and the invertible integral H{z) = 
J [b(z)]~ 1 dz exists on the domain of definition of the original equation (|7.1|l for X. 
In this case one finds 

X 1 {t) = H-\W{t)+H{X 1 {Q)))- 

Once an appropriate splitting is found then it is relatively easy to implement the 
related numerical algorithm. The proposed splitting-step method efficiently works 
since its implementation essentially relies on the well-known variation-of-parameters 
formula for perturbed dynamical systems which extends to SDEs. Recall that, by this 
formula, if the equation 

dX(t) = /3(X (t),t)dt + a(X (t),t)dW(t) 

has known fundamental solution <!> = &(t,Xo) then the exact solution of the original 
equation l|2.1|l possesses the pathwise representation 

/t+At 
[$(s,X(t))]- 1 a(X(s),s) ds 

w $(f + At, X(t)) + a(X(t), t)At 

on each subintervals [t, t + At] C [0,T]. Thus, the motivation of our splitting- 
step technique is apparent by finding $ and numerical integration of expressions 
/ a(X(s), s) ds. 
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Appendix A. Non-central chi-square distribution random number gen- 
eration. The splitting-step method proposed in this paper relies on the exact nu- 
merical sampling of transition probability density for some processes. In particular, 
efficient generation of random numbers x'd W with a non-central chi-square distribu- 
tion with d degrees of freedom whose probability density function is found in ^B] by 
noting that 

-(A+x)/2 . . (d-2)/4 , 

(A.l)V[ x , d 2 (\)=x]=p(x;d,\) = — (-) / ( „-a)/2(VA^), x > 0. 

This distribution is properly defined for any d positive, and was extended to the case 
d = by Siegel [2S] • Here we will extend it to the case d = — 2, — 4, . . . and will show 
how to sample this distribution. 

To this end, we use the fact that the distribution l|A.l|) can be expressed also as 
a mixture of central x 2 variables with Poisson weights 

(A.2) p( x;d; A)=]T^^ p (x;d + 2j), d > 0, 

where po(x; d) is the distribution of a chi-square random variable x\ with d degrees of 
freedom. This expression suggests a simple and efficient procedure to obtain XdW 
random variables: 

1. Choose K from a Poisson distribution with mean A/2 so that V[K = k] = 
e- k ' 2 {\/2) k /k\ (A: = 0,1,...). 

2. Then take X^ 2 (A) = Xd+2K> which can be done using the any standard random 
number generator of the x\ distribution. 

In the d = case, the x'dW distribution has a discrete component at zero with 
mass e A / 2 (which represents the probability to get absorbed at zero in our stochastic 
processes), see [35]. We have 

(A.3) p (a; ;0,A) = ]ri^ p (x;2j)+e- x / 2 6(x), d = Q. 

i=i J ' 

The procedure above can be modified to account for this discrete component by taking 
the convention that the central x\ distribution is identically zero when d = 0. This 
convention can be extended to even negative dimensions to get 

(A.4) p{x;d,\)= (A/2)J r V2 Po(z; d + 2j) + S(x) 'f) (W^ 2 

] = \d\/2+l J ' j=0 3 ' 

for d = 0, -2, -4, . . .. 

Summarizing, if K is a Poisson random number with mean A/2 we have 

(A.5) XdW=Xd+2K, d>0 

and 

( ) X " (A) = \xW ifd + 2if;0 ^ = 0,-2,-4,.... 
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This sampling of the probability distribution function is exact and should be used 
especially when A is small. However, when A is large, the XdW distribution asymp- 
totically converges to the Gaussian distribution and other approximations (like the 
ones in ^5]) might be considered to improve the speed of our algorithm. 

REFERENCES 

[1] M. Broadie and O. Kaya, Exact Simulation of Stochastic Volatility and other Affine Jump 

Diffusion Processes, Working Paper (2004). 
[2] J. C. Cox, Notes on option pricing I: Constant elasticity of variance diffusions, Working Paper, 

Standford University (1975); reprinted in Journal of Portfolio Management 22 (1996), pp. 

15-17. 

[3] J. C. Cox, J. E. Ingersoll, and S. A. Ross, A Theory of term structure of interest rates, Econo- 

mctrica 53 (1985), p. 385-407. 
[4] F. Castell and J. Gaines: The ordinary differential equation approach to asymptotically efficient 

schemes for solution of stochastic differential equations, Ann. Inst. H. Poincare Probab. 

Statist. 32 (1996), p. 231-250. 
[5] D. A. Dawson and E. A. Perkins Measured-valued processes and renormalization of branching 

particle systems, in Stochastic Partial Differential Equations: Six Perspectives, R. A. Car- 

mona, B. Rozovskii eds. (AMS Mathematical Surveys and Monographs, No. 64., 1998), p. 

45-102. 

[6] H. Doss: Liens entre equations differentielles stochastiques et ordinaires, Ann. Inst. Henri 

Poincare XIII (1977), Section B, No. 2, p. 99-125. 
[7] R. Durrett and E. A. Perkins: Rescaled contact processes converge to super- Brownian motion 

in two or more dimensions, Probab. Theory Relat. Field 114 (1999), p. 309-399. 
[8] E. B. Dynkin, Superprocesses and partial differential equations, Ann. Probab. 21 (1993), p. 

1185-1262. 

[9] W. Feller, The general diffusion operator and positivity preserving semi-groups in one dimen- 
sion, Ann. Math. 60 (1954), p. 417-436; Diffusion processes in one dimension, Trans. 
Amer. Math. Soc. 77 (1954), p. 1-31; The parabolic differential equations and the associ- 
ated semi-groups of transformations, Ann. Math. 55 (1952), p. 468-519. 

[10] W. Feller: Probability theory and its applications, 2nd edition, Wiley, New York, 1971. 

[11] J.G. Gaines, Numerical experiments with S(P)DE's, in Stochastic Partial Differential Equa- 
tions, A.M. Etheridge. ed., London Mathematical Society Lecture Note Series 216, Cam- 
bridge Univ. Press., 1995, p. 55-71. 

[12] D.J. Higham and X. Mao: Convergence of Monte Carlo simulations involving the mean- 
reverting square root process, Journal of Computational Finance 8 (2005), p. 35-62. 

[13] D.J. Higham, X. Mao, and A.M. Stuart: Strong convergence of Euler-type methods for nonlin- 
ear stochastic differential equations, SIAM J. Numer. Anal. 40 (2002), p. 1041-1063. 

[14] H. Hinrichsen N on- equilibrium critical phenomena and phase transitions into absorbing states, 
Advances in Physics 49 (2000), pp. 815958. 

[15] N. L. Johnson, S. Kotz, and N. Balakrishnan, Continuous Univariate Distributions (Wiley, New 
York, 1994), Vol. II; 

[16] C. Kahl and H. Schurz: Balanced Milstein methods for ordinary SDEs, Preprint M-05-004, 

SIU, Department of Mathematics, Carbondale, 2005. 
[17] P.E. Kloeden and E. Platen: Numerical solution of stochastic differential equations, Springer, 

New York, 1998. 

[18] P.E. Kloeden, E. Platen, and H. Schurz: Numerical solution of SDEs through computer exper- 
iments, Springer, New York, 1994. 

[19] G.N. Milstein: Numerical integration of stochastic differential equations, Kluwer, Dordrecht, 
1995. 

[20] G.N. Milstein, E. Platen, and H. Schurz: Balanced implicit methods for stiff stochastic systems, 
SIAM J. Numer. Anal. 35 (1998), p. 1010-1019. 

[21] T. Misawa: A Lie Algebraic approach to numerical integration of stochastic differential equa- 
tions, SIAM J. Sci. Comput. 23 (2001), p. 866-890. 

[22] E. Moro: Numerical schemes for continuum models of reaction- diffusion systems subject to 
internal noise, Phys. Rev. E 70 (2004), p. 045102-1,045102-4. 

[23] C. Mueller and R. Tribe: A phase transition for a stochastic PDE related to the contact process, 
Probab. Theory Relat. Fields 100 (1994), p. 131-156; Stochastic p.d.e.'s arising from the 
long range contact and long range voter process, Probab. Theory Relat. Fields 102 (1995), 



NON-NEGATIVITY PRESERVING NUM. ALGORITHMS FOR SDEs 



23 



p. 519-545. 

[24] L. Pechenik and H. Levine, Interfacial velocity correction due to multiplicative noise. Phys. 
Rev. E 59 (1999), p. 3893. 

[25] H. Schurz: Numerical regularization for SDEs: Construction of nonnegative solutions, Dynam. 
Syst. Appl. 5 (1996), p. 323-352. 

[26] H. Schurz: Stability, stationarity, and boundedncss of some implicit numerical methods for 
stochastic differential equations and applications, Logos- Verlag, Berlin, 1997. 

[27] H. Schurz: Numerical analysis of SDEs without tears (invited chapter), in Handbook of Stochas- 
tic Analysis and Applications, V. Lakshmikantham and D. Kannan (eds.), Marcel Dekker, 
Basel, p. 237-359, 2002. 

[28] H. Schurz: General theorems for numerical approximation of stochastic processes on the Hilbert 

space H 2 ([0,T]), Electr. Trans. Numcr. Anal. 16 (2003), p. 50-69. 
[29] A. F. Siegel, Biometrika 66 (1979), p. 381. 

[30] G. Slade: Scaling limits and Super- Brownian Motion, Notices to the AMS 49 (2002), p. 1056. 

[31] W. Wagner and E. Platen: Approximation of Ito integral equations, Technical Report at ZIMM, 
Acad. Sci. GDR, Berlin, February, 1978 (obtainable from library at WIAS, Berlin, Ger- 
many, http://www.wias-berlin.de/ I. 



